D = zeros(dimU+dimV,dimU+dimV);
for i=1: NIut
    for j=1:NJut
        if (i<=NIue && j>NJue) continue; end
        
        indexU_cen   = Umap (i  ,j  );
        
        if (j<=NJue)
            
            % Boundary conditions on inlet & outlet
            if (i==1) continue ; end
            if (i==NIut) continue ; end    % Boundary condition oulet - null normal gradient
            
        else
            
            % Boundary conditions on expansion & outlet
            if (i==NIue+1) continue ; end    % Boundary condition expansion
            if (i==NIut) continue ; end    % Boundary condition oulet - null normal gradient
            
        end
        
        
        %% Dx coefficients
        centerCoeff = 4;
        indexU_cen   = Umap (i  ,j  );
        
        indexU_right = Umap (i+1,j  );
        D(indexU_cen,indexU_right) = 1/h^2;
        indexU_left  = Umap (i-1,j  );
        D(indexU_cen,indexU_left) = 1/h^2;
        
        
        
        
        % Top and bot volumes, considering ghost nodes... The ghost nodes
        % they would add, isn't it?? because they are -center
        if ( j==1 )
            centerCoeff = centerCoeff + 1;
        else
            indexU_bot   = Umap (i  ,j-1);
            D(indexU_cen,indexU_bot) = 1/h^2;
        end
        if (j<=NJue)
            if ( j==NJue && i<=NIue )
                centerCoeff = centerCoeff + 1;
            else
                indexU_top   = Umap (i  ,j+1);
                D(indexU_cen,indexU_top) = 1/h^2;
            end
        else
            if ( j==NJut )
                centerCoeff = centerCoeff + 1;
            else
                indexU_top   = Umap (i  ,j+1);
                D(indexU_cen,indexU_top) = 1/h^2;
            end
        end
        
        D(indexU_cen,indexU_cen) = -centerCoeff / h^2;
    end
end


for i=1:NIvt
    for j=1:NJvt
        if (i<=NIve && j>NJve) continue; end
        %% Dy coefficients
        centerCoeff = 4;
        indexV_cen  = dimU + Vmap (i  ,j  ) ;
        
        if (j<=NJve)
            % Boundary conditions on top and bot boundaries
            if (j==1)  continue ; end
            if (j==NJve && i<=NIve ) continue ; end
        else
            % Boundary conditions on top and bot boundaries
            if (j==NJvt )  continue ; end
        end
        
        if (i==NIvt)
            centerCoeff = centerCoeff ;%+ 1; ---- > convective outflow modification
        else
            indexV_right = dimU + Vmap (i+1,j );
            D(indexV_cen,indexV_right) = 1 / h^2;
        end
        if (j<=NJve)
            % Right and left boundaries, considering ghost nodes
            if (i==1)
                centerCoeff = centerCoeff + 1;
            else
                indexV_left = dimU + Vmap (i-1,j );
                D(indexV_cen,indexV_left) = 1 / h^2;
            end
        else
            % Right and left boundaries, considering ghost nodes
            if (i==NIve+1)
                centerCoeff = centerCoeff + 1;
            else
                indexV_left  = dimU + Vmap (i-1,j );
                D(indexV_cen,indexV_left) = 1 / h^2;
            end
        end
        indexV_bot =  dimU + Vmap (i,j-1);
        D(indexV_cen,indexV_bot) = 1 / h^2;
        indexV_top =  dimU + Vmap (i,j+1);
        D(indexV_cen,indexV_top) = 1 / h^2;
        
        D(indexV_cen,indexV_cen) = - centerCoeff / h^2;
        
    end
end
